Surrogate-based optimization design for surface texture of helical pair in helical hydraulic rotary actuator

A good surface texture design can effectively improve the tribological performance of the helical pair within a helical hydraulic rotary actuator(HHRA). However, the optimization design process can be time-consuming due to the multiple design variables involved and the complexity of the mathematical model. This paper proposes a modified efficient global optimization (MEGO) method for solving such demanding surface texture design challenges. The MEGO utilizes a Kriging model with the optimized Latin hypercube sampling (OLHS) for initial sampling and the proposed modified expected improvement (MEI) function for sequential sampling. A comparative study of several global optimization algorithms with the MEGO on the surface texture design is performed. Subsequently, surrogate-based optimization and parameter analysis are carried out, resulting in the identification of an optimal set of texture parameters. The findings reveal the superiority of the MEGO in both model prediction accuracy and refinement of minima. Moreover, compared to the base design, the friction coefficient can be reduced by up to 45.2%.

helical pair is subjected to mixed lubrication, necessitating the consideration of asperity contact.The helical shape of the oil film in a helical pair further necessitates transformation into an equivalent plane, introducing additional structural parameters and increasing the complexity of the mathematical model.Lastly, the oil film thickness between helical pairs, primarily dependent on the external load, is yet to be definitively determined and requires numerous iterative calculations, rendering the design procedure more time-consuming and challenging to execute.
With their independence from gradient information and ease of implementation, meta-heuristic algorithms have emerged as an efficacious approach to texture optimization.Wang et al. 3 used a hybrid evolutionary algorithm for the optimization design of texture bottom contours, demonstrating that global-optimum bottom contours exhibited superior load capacity compared to regular ones.Similarly, Zhang et al. 4 introduced an optimal design scheme for surface texture based on a genetic algorithm (GA), successfully enhancing the unidirectional sliding tribological performance.Their results indicated that fish and bullet shapes delivered optimal tribological performance.Chen et al. 5 utilized a nondominated sorting genetic algorithm II (NSGA-II) to optimize the surface texture at the interface between the cylinder and valve plate in an axial piston pump, with simulation results showing that water-drop-shaped dimples displayed superior tribological behaviors and could reduce leakage by an average of 9.9%.Tang et al. 6 employed a multi-objective hybrid evolutionary algorithm to maximize load carrying capacity while minimizing the friction coefficient in an axial piston pump by optimizing the surface texture of the slipper bearing, achieving improvements of 64.8% and 4.5% respectively after optimization.Bei et al. 7 utilized the average Reynolds equation to solve the rough contact model and employed GA to ascertain the optimal texture shape.While these methods have been widely and successfully implemented, they necessitate extensive calculations and iterations during the optimization process.This drawback renders them suitable only for the optimization of surface texture that involves short-duration simulations and a limited number of design variables.
Surrogate modeling methods, such as support vector regressions (SVR) 8 , Kriging models 9,10 , and radial basis functions (RBF) 11 , employ fast-computation surrogate models to fit time-consuming simulations, thereby significantly reducing the optimization process.Within the realm of engineering, there exist many successful instances where the approximate method is employed to simulate expensive numerical simulations for model prediction [12][13][14] , subsequently leading to the proliferation of numerous surrogate-based optimization algorithms [15][16][17] .Among them, the efficient global optimization (EGO) algorithm 18 is frequently utilized due to its high efficiency and global optimization capability with limited sample points.The EGO is one of the most widely studied and applied surrogate-based optimization algorithms, with its applications primarily in aerodynamics.Jeong, Murayama, and Yamamoto 19 introduced a novel global optimization method that combines EGO and genetic algorithm, suitable for aerodynamic design.Ghoreyshi, Badcock, and Woodgate 20 presented a hybrid sampling strategy that combines maximizing the mean square error(MAXMSE) criterion and expected improvement (EI) function within the EGO framework, which was applied to the generation of aerodynamic tables for flight simulation.Ariyarit et al. 21demonstrated a hybrid surrogate-based optimization method premised on standard EGO for aerodynamic optimization, showcasing the capacity of the method to generate superior blade shapes with improved aerodynamic efficiency.Guzman Nieto, ElSayed, and Walch 14 introduced an EGO-based algorithm to effectively predict critical dynamic aeroelastic loads, with the results indicating a total time reduction of 40.79% compared to the sole use of Kriging interpolation.Despite EGO's commendable performance in solving these optimization problems, its deployment in texture design necessitates further exhaustive investigation.Certain constraints limit its usefulness: EGO primarily centers around local search, potentially leading to convergence to a local optimum as opposed to a global one when encountering potent nonlinear problems 22 .This focus can also result in poor model precision when sample points are limited in high-cost computation scenarios.Thus, the development and application of a suitable surrogate-based optimization algorithm for the multi-variable, strongly nonlinear model of surface texture warrants further exploration.
To reduce the friction of helical pair in HHRA by reasonable surface texture design, and improve the efficiency of the optimization design process, a modified efficient global optimization method, which takes into account the shortcomings of the standard EGO algorithm, is proposed.The first step in this method constitutes the construction of an initial Kriging model, with the OLHS strategy being employed to procure a set of uniform sample points for the establishment of the Kriging model.The second step involves the addition of new infilling samples to seek the global minimum and enhance model precision with minimal expenditure; a modified expected-improvement sampling criterion, capable of adaptively selecting the infilling samples in accordance with the optimization process, is utilized as the sequential sampling approach.Following this, the mathematical model of the textured helical pair under the mixed lubrication regime is established, based on which the proposed MEGO is utilized to construct a surrogate model for the low friction coefficient design of the surface texture.To demonstrate the superiority of MEGO in surface texturing design, a comparative study involving several global optimization algorithms and the MEGO with respect to surface texture design is executed.Lastly, surrogate-based optimization and analysis are conducted, resulting in the identification of an optimum set of texture parameters of the elliptical-shape dimples.Parameter analysis is then undertaken using the global sensitivity analysis method, correlation analysis, and controlled variable method.

Optimized Latin hypercube sampling
In 1979, McKay, Beckman, and Conover 23 originally proposed Latin hypercube sampling(LHS).Afterward, Stein 24 proposed a mathematical approach to this method, it can be expressed as an N × M matrix Q: where N is the number of selected points, Q i is the ith sample point, q ij is the value of the jth dimension of the ith sample point.
LHS employs multi-dimensional stratified sampling to fill the entire design space without overlap.Despite this, LHS suffers from inadequate spatial uniformity, making it challenging for initial sample points to capture the global information of the actual model.To address this, Liefvendahl and Stocki 25 devised an optimization criterion function to enhance the distribution uniformity of samples.The objective function of this optimized Latin hypercube sampling (OLHS) criterion takes the following form:

Kriging model
The Kriging model, a renowned surrogate model based on the Gaussian process, has found widespread use across numerous engineering fields, particularly following its application to computer simulation by Sacks et al. 26 .An ordinary Kriging model can be articulated as follows: where µ denotes the mean value of the Kriging prediction function, while ǫ(x) is the error term of a Gaussian process with zero mean and nonzero variance.The covariance of this term can be defined as: where c x i , x j is a correlation function primarily influenced by the distance between the sample points x i and x j , and R() denotes the symmetric correlation matrix.
While various forms can be adopted for φ x i , x j , including exponential, Gaussian, cubic, and others, this study employs the most commonly used Gaussian form: where num signifies the number of variables, and θ k represents a correlation parameter, previously unascertained, that can be computed by resolving the maximum likelihood function.
The predictive value f for untried sample point x within the Kriging model can be conveyed as: The Kriging variance, also known as the mean squared error (MSE), is expressed as follows: (1)

Modified expected improvement function
The EI function was initially introduced as a statistical criterion for Kriging-based optimization by Jones, Schonlau, and Welch 18 .This function leverages the Kriging model to glean statistical information from each sample point, guiding the selection of new infilling points and the overall search direction of the EGO algorithm.
The EI function can be expressed as follows: where f min signifies the minimum value of the real function.ψ(•) represents a probability density function, while ξ(•) symbolizes a cumulative distribution function with a standard normal distribution.The expression ( f min − f (x) indicates a potential improvement, facilitating a more localized exploitation of the search space, whereas s(x) reflects the uncertainty inherent in the surrogate model, which is responsible for a more globally explorative search, quite the opposite.Consequently, the EI function can seamlessly balance global exploration and local exploitation within the optimization process.While the EI function can locate an optimal value in certain simplistic cases, the computed standard deviation s(x) tends to be marginally lower than the actual value, a phenomenon often referred to as "underestimation" 22 .Consequently, only points in close proximity to the current optimal solution register significant EI values, leading to a more detailed exploration of the adjacent region until the local uncertainty decreases sufficiently, prompting a more global search.
One way to utilize the EGO framework for a more accurate surrogate model and more global search in the early stage is to modify the EI function.Sóbester, Leary, and Keane 27 presented a weighted expected improvement criterion allowing for a more flexible means of balancing global and local search results.Ponweiser, Wagner, and Vincze 22 presented the "generalized expected improvement" (GEI) criterion, founded on the EI criterion, where the parameter "g" determined whether the algorithm tends to local exploitation or global exploration.Despite the enhanced global search capabilities of the EGO algorithm offered by these methods, they still require manual parameter setting.In practical applications, selecting an optimal value for sampling poses a considerable challenge for designers.
Given the slight underestimation of the prediction error by the EGO algorithm in Formula9, the global reach of the EI function can be extended by incorporating a weighting factor preceding the term s(x) .This factor should evolve with the iterations, as a global exploration at the onset of the EGO algorithm is preferred, while local exploitation is desirable in the latter half 22 .
To fulfil this requirement, the smooth, stepwise increasing behaviour of the sigmoid function is leveraged, and a variant sigmoid function 28 is introduced and rewritten as the trend term for the weight coefficient.This modification enhances the global scope in the first half of the iteration, and allows the weight coefficient to gradually decrease with the iteration of the sample point.The formula is expressed as follows: where n iter signifies the total number of sample updates, whereas i iter represents the i iter th update.
The correlation between w iter and i iter under varying values of n iter is depicted in Fig. 2. It can be observed that w iter approximates 1 in the initial third of n iter and swiftly transitions from 1 to 0 in the regions spanning from the first third to half of n iter .In essence, the global weight coefficient escalates in the first half of n iter , facilitating a more global search.For the latter half of w iter , the weight value is nearly 0, implying that the EI algorithm is retained, thus ensuring a more localized search in the second half of the optimization process.
Additionally, θ k in Formula5 emerges as a critical parameter of the Kriging model.A larger θ k value augments the fluctuation of the Kriging model curve, while a smaller θ k value smooths the curve.Greater values of θ k sug- gest the need for a more global search, and the converse is also true.Therefore, θ k serves as an excellent choice for a weighting regulator to adjust the weight ratio.As a vector, θ k can be quantified via the introduction of the Euclidean norm and normalized using the arctan function.
Applying these principles, the weighting regulator w θ k can be formulated as follows: Therefore, the MEI function can written as: www.nature.com/scientificreports/ The MEI is correlated with i iter and θ k , changing in accordance with updates to the Kriging model.This correlation facilitates a more global search during the initial stages and a more localized search near the end of an iteration, with the weight ratio adapting in tandem with updates to the Kriging model.

Geometry and governing equations
The equivalent section of the interface of the helical pair, as depicted in Fig. 3, is derived using the method proposed by El-Sayed and Khatan 29 .The geometric relationship between the actual model and the equivalent model can be expressed as follows: where symbol ˜ represents the variables are located on the equivalent plane, θm is the equivalent angle from the start to the end of helicoid.r i is the inside radius of nut, r o is the outside radius of screw, and Figure 4 shows the cross-section of the screw tooth.Based on the equivalent model, the axial film thicknesses at any radius r can be written as follows: where h a0 is the assumed axial clearance, h is the clearance of screw and nut along the axial direction, the sub- scripts t/b represents the oil film on the top/bottom side.r = arctan(P/(2πr)) , where r denotes the helix angle at any radius r, P denotes the screw lead.α r = arctan (tan α cos r ) , in which α is pressure angle.
Building on the findings of Yu, Wang, and Zhou 30 , elliptical-shaped dimples deliver the optimal load-carrying capacity compared to other regular-shaped dimples such as circles, triangles, and rectangles.As a result, this study employs elliptical-shaped dimples as the textured structure of the helical pair.The representative section of the interface of the helical pair with elliptical-shaped dimples is illustrated in Fig. 5, where the dimples are uniformly distributed in both radial and circumferential directions.
To facilitate calculation, it is assumed that the texture is positioned within a rectangular unit of dimensions l × w .The parameters defining the dimples are as follows: Z c represents the number of circumferential dimples, Z r the number of radial dimples, and φ the orientation of the major axis of ellipse.B denotes the length of the texture area in the radial direction, while θ t signifies the length of the texture area in the circumferential direc- tion.The dimple depth is h d , with l a and l b specifying the major and minor axes of the ellipse, respectively.The spacing between dimples in the radial (l) and circumferential (w) directions, the area ratio S p , the axial ratio of the ellipse r t , and the texture region can be respectively expressed as follows: Evidence has shown that provided the Reynolds number is adequately low and the aspect ratio of the dimples is sufficiently large, the Reynolds equation proves effective 31 , a condition applicable to our model.The assumption is that the oil film represents an incompressible Newtonian and laminar flow.Consequently, the equivalent Reynolds equation 32 can be expressed as:  where � r = 1.5 − 5 sin 2 r + 5 sin 4 r / 60 cos 2 r , in this context, r and θ are the coordinate directions of the polar coordinates õ − r θ , situated on the equivalent plane of the helicoid.The hydrodynamic pressure is denoted by p, with ρ representing the oil density and η its viscosity.The rotational speed is given by ω , and z nut signifies the axial displacement disturbances.The subscript t/b refers to the top side/bottom side.
Generally, the HHRA operates under high load and low-speed conditions, resulting in a relatively thin film between the helical pair, necessitating the consideration of asperity contact force.This paper employs the Greenwood and Tripp model 33 , which is widely used in asperity contact studies.By assuming a Gaussian distribution for the height of the micro-convex bodies on the friction surface, the asperity contact pressure can be defined as: where , where E 1 and E 2 denote the elastic modulus of the two bodies in contact, and E ′ is the composite elastic modulus.The Poisson's ratios of the contacting bodies are represented by v 1 and v 2 .The asperity curvature radius is indicated by β , while A signifies the asperity density.The composite surface roughness is designated by σ , and H σ is defined as the ratio of the mean film thickness to σ.
The statistical function F n (u) can be written as: where s is the correlation parameter, set at 6.804 33 .
Implementing the average flow model 34 , the equivalent Reynolds equation that takes into account roughness effects can be articulated as: where φ x and φ y represent the pressure flow factors along the x and y directions, as defined by Patir and Cheng 34,35 .The contact factor is denoted by φ c , while φ s signifies the shear flow factor, as outlined by Wu and Zheng 36 .The lead of the screw/nut is represented by P.

Velocity boundary
Utilizing velocity decomposition, the velocities along the top and bottom sides of the screw tooth can be determined as follows: where v θ1/ θ2 and v z1/z2 denote the speeds in the θ and z directions, respectively.

Pressure boundary
The external pressure of the helical pair equals 0.1 MPa, and the pressure between the contact surface of the helical pair corresponds to the oil film load capacity.Besides, the minimum value of cavitation pressure is specified as 0.03 MPa 5 .

Cavitation boundary conditions
Owing to its high accuracy and ease of implementation, the Reynolds cavitation boundary condition has been effectively utilized in numerous studies investigating surface texture 4,5,30,37 .Therefore, the Reynolds cavitation boundary condition is chosen in this paper.

Calculation of friction coefficient
The expression of friction coefficient can be written as: where, the total load capacity D is an amalgamation of the load-carrying capacity induced by hydrodynamic support and that induced by asperity contact.This can be calculated as follows: The total friction force f comprises both the hydrodynamic friction force and the asperity contact friction force.It can be represented as: where v denotes the linear velocity, while φ f and φ fs indicate the flow factors induced by friction, as defined in the studies by Patir and Cheng 34,35 .The asperity contact friction coefficient is represented by µ asp .
For the HHRA studied in this paper, the oil film thickness between the helical pair remains undetermined.Its value is primarily influenced by the external load driven by the HHRA, necessitating computation from the external load.The comprehensive calculation procedure is depicted in the "oil film characteristics module" in Fig. 6.
The parameters used in the numerical analysis are shown in Tables 1 and 2.
The equivalent Reynolds equation that incorporates roughness effects is discretized utilizing the finite difference method (FDM).To calculate the film pressure distribution, the successive over relaxation (SOR) Gauss-Seidel iterative method, with an over-relaxation factor at each discretized node, is employed.The over-relaxation (20)     factor is set between 1.5 and 1.8 37 , with a chosen value of 1.6 in this study to expedite convergence.MATLAB R2022b software facilitates the construction and resolution of the equivalent Reynolds equation, operating on a computer CPU that is i7-9750h with 16 G RAM. Grid density significantly influences the accuracy of results.While fine-meshed grids can enhance computational accuracy, they also demand more computing time.Table 3 presents the computation time of this mathematical model at varying grid densities.Relative to the finest grid, a grid density of 120 × 720 is selected for texture optimization, as it offers the shortest computing time when the error is within 1%.However, as the computation with the 120 × 720 grid remains time-intensive, it is employed solely in the final optimization(5.2) and sensitivity analysis (5.3).For the comparison of algorithms (5.1), which necessitates the generation of a considerable number of sample points, the 40 × 240 grid is chosen.This selection is due to a greater emphasis on the alignment of the predicted value with the actual one, and the ability to locate the minimum value, rather than obtaining an accurate friction coefficient and optimal texture parameters.

Procedures of modified efficient global optimization
A flowchart detailing the optimization design process for the surface texture of the helical pair is provided in Fig. 6.A comprehensive explanation follows.
Step 1: Establish the optimization objectives, the ranges for design variables, and the objective and constraint functions.
Numerous studies indicate that within a certain range, the friction coefficient progressively decreases with the increasing textured area density.However, exceeding this range could elevate contact stress, limiting the reduction of the friction coefficient and intensifying wear.Xu et al. 38 found that elliptical-shaped dimples with an area density of 10.6-14.1% yielded the best performance.To minimize the friction coefficient, a 14% area density was selected for this study.The remaining undefined geometry and distribution parameters of the elliptical-shaped dimples serve as the design parameters.Notably, most surface texture optimization research employs dimple depth and diameter as design variables for circular dimples, given their significant impact on the lubrication characteristics of surface texture.For a single texture unit, the area density directly correlates with the square of the dimple diameter.Thus, the dimple depth and diameter as design variables sufficiently encompass all necessary design information.Yet, the interaction among the dimples can also influence hydrodynamic lubrication for multiple texture units.Consequently, the axial number and the circumferential number of the dimples are chosen as part of the design variables.When the area density remains constant, these numbers not only determine the dimple diameter but also reflect the spacing between the dimples in the radial and circumferential direction.
The remaining undefined geometry and distribution parameters of the elliptical-shaped dimples function as the design parameters, with their range chosen to be technologically feasible and prevent dimple overlap.
The optimization issue concerning the surface texture of the helical pair, as considered in this study, can be expressed mathematically as follows: www.nature.com/scientificreports/ Step 2: Conduct a design space exploration and sample initial points using the design of experiment (DOE) method.In this study, the OLHS method mentioned above is used to obtain maximum information by uniformly adding sample points from the design space.
The initial number of sample points is determined by the number of variables.According to Nuchitprasittichai and Cremaschi 39 , the initial number of sample points is typically ten times the number of variables.Considering that this study involves five design variables, an initial sample count of 50 is chosen.However, due to the constraints of the OLHS method, sample points are not selected at the boundary, resulting in an incomplete coverage of the design area.The boundary points are defined by upper and lower bounds, and with five variables included in this study, the total number of boundary points is 2 5 .As a result, an aggregate of 82 (50 + 2 5 ) initial sample points is generated.
Step 3: Conduct the numerical analysis as outlined in "Mathematical model", using the initial sample points determined in Step 2.
Step 4: Use the results from Step 3 to establish the initial Kriging model, representing the mapping relationships between the objective functions and the design variables.
Step 5: Execute the sequential sampling strategy, employing the MEI function to boost the accuracy of the Kriging model and locate the minimal values of the actual model through the addition of new sample points.The sequential sampling process concludes when the count of iteration steps for sequential sampling attains the predetermined value.
To evaluate model accuracy, 50 sample points are used to calculate the root mean square error(RMSE) value of the Kriging model.The stopping criterion during the sampling process is the total number of sample points reaching 200 times.Therefore, 118 iterations are required, given that the initial number of sample points is 82.
Step 6: The final Kriging model, constructed using both initial and sequential sampling data, serves to predict the targets, effectively replacing the need for laborious numerical analysis, and enables the procurement of optimal design results.

Evaluation index of surrogate-based optimization algorithm
Predictive accuracy and optimization potential are critical factors for surrogate-based optimization algorithms.To appraise the accuracy of the surrogate model in quantitative terms, the introduction of three evaluation indices is proposed as follows: www.nature.com/scientificreports/where f predict denotes the predicted value obtained from the surrogate model, while f actual refers to the actual value derived from numerical analysis.The optimization potential is evaluated as follows: where opt.predict refers to the globally optimized solution obtained from the surrogate-based optimization algorithm, while opt.actual denotes the actual global optimized solution.

Algorithm comparison
To demonstrate the superiority of the MEGO method in surface texture design, it is compared with commonly used surrogate-based optimization algorithms: EGO 18 , Kriging+GA 16 , RBF+GA 15 , BP+GA 40 , and the typical meta-heuristic algorithm GA.
It is worth noting that the algorithm comparison, distinct from the optimization steps outlined in 4.1, prescribes a fewer number of initial sample points and a maximum limit for sample points.This arrangement is designed to evaluate their ability to seek optimization and fitting accuracy under the restriction of limited sample points.Hence, each global optimization algorithm is limited to 60 points, equivalent to 60 simulation times.Initially, 20 points are generated using the OLHS method, and the remaining 40 infilling points are selected by the MEI function for each iteration.To test the accuracy of the established surrogate model, 100 random test samples are generated.Given the randomness of the sampling distribution induced by DOE, each global optimization algorithm is computed 20 times for every case, and the results are averaged.The standard deviation (std) is employed here to gauge the stability of the optimization results.For ensuring the accuracy of the prediction model, MRE and MaxRE should be as minimal as possible, and a lower RMSE evaluation index signifies a more accurate prediction model.Figure 7 displays the test points predicted by each algorithm, and the comparative ( 25) www.nature.com/scientificreports/results are depicted in Table 4 and Fig. 8. Notably, GA, while capable of global optimization, lacks model prediction ability, hence it is excluded from discussions about the accuracy of the prediction model.The minimum µ cof is 0.0433, as calculated by the Monte Carlo sampling method.It is clearly observed that the Kriging-based optimization algorithms (MEGO, EGO, Kriging + GA) showcase lower MRE, MaxRE, and RMSE values compared to the others (RBF + GA, BP + GA).This highlights the superior ability of the Kriging-based surrogate model to fit the nonlinear curve of the actual model in this case.Moreover, the MRE, MaxRE, and RMSE values for EGO are marginally higher than those for MEGO, signaling the superior prediction ability of MEGO over the standard EGO.
Regarding the optimal value, the global optimization ability is sorted according to the magnitude of (µ cof ) min as follows: MEGO > EGO > Kriging + GA > GA > BP + GA > RBF + GA.It is evident that the global optimization ability of MEGO surpasses that of the other models, being only 1.4% away from the optimal result.In contrast, EGO, Kriging + GA, RBF + GA, BP + GA, and GA are 3.0%, 15.7%, 69.1%, 43.4%, 24.2% away from the optimum result, respectively.This implies that MEGO can yield the optimal global solution with the fewest iterations.Furthermore, MEGO exhibits the smallest standard deviation among all global optimizations, indicating that MEGO consistently delivers the most stable optimization results.

Optimum texture parameters
Before embarking on the parameter analysis, it is crucial to confirm that the accuracy demands are fulfilled by the constituted Kriging model.At the same time, it must be assured that a local minimum does not stand as the final optimization outcome in the course of texture design optimization.To support this, the RMSE values from each step of the sampling process, obtained via the MEI method, together with the minimum friction coefficient values, are recorded in this paper.A complete visual representation is available in Fig. 9a,b.
For the friction coefficient studied in this paper, its value ranges from 0.04 to 0.09.When the RMSE value is less than 0.003, the established surrogate model can be considered to meet the engineering accuracy requirements.As discernible from Fig. 9a, the RMSE values of the Kriging model at a rotational speed less than 60 rpm remain below 0.003(the RMSE value is still less than 0.003) at all sample points, exhibiting no significant fluctuations with the increase in iteration steps (the RMSE value is still less than 0.003).This suggests that a highly accurate Kriging model can already be established from the initial sample points, with no overfitting observed upon the addition of more sample points.Moreover, the RMSE values of the Kriging models at rotational speeds of 90rpm and 120rpm initially exhibit considerable fluctuations and are greater than 0.003.However, as the number of iteration steps increases, the RMSE values of both models gradually converge to less than 0.003 after about 175 sample points and remain so at the 200th step.Hence, it can be concluded that when the number of The optimal texture parameters are determined using the MEGO.The base design, built on experiential knowledge, acts as a reference group.The optimal design, reference group, smooth surface, and corresponding friction coefficients at various rotational speeds are showcased in Table 5.The percentage improvement in the friction coefficient after optimization is illustrated in Fig. 10.The pressure distributions and texture profiles of the optimal design and its reference group are depicted in Figs.11 and 12, respectively.
As indicated in Fig. 11a, rapid shifts in pressure distribution can be observed in the vicinity of the dimples.The peak pressure is found at the starting position of the dimples along the sliding direction, a result of the additional hydrodynamic pressure from the converging wedge within this region.Conversely, a diverging wedge appears at the terminal position of the dimples along the sliding direction, inducing a negative pressure where cavitation transpires due to the Reynolds boundary conditions.Given the larger area of the positive pressure zone compared to the negative pressure zone and the influence of Reynolds boundary conditions, the textured area can provide an additional hydrodynamic load-bearing capacity.However, as shown in Fig. 12a, the additional hydrodynamic effects are not pronounced, and load-bearing is primarily facilitated by asperity contact.This leads to a reduction in film thickness and subsequently a gradient descent of pressure distribution along the radial direction across the entire surface, attributable to varying helix angles at the corresponding radius.
Table 5 demonstrates that all optimal texture parameters contribute positively to the reduction of the friction coefficient of the helical pair, with the extent of influence being associated with rotational speed.This can be attributed to the significant impact of sliding movement on the aforementioned additional hydrodynamic load-bearing capacity, as displayed in Fig. 11a.A certain range of high rotational speed can enhance this capacity.Moreover, the friction coefficient of the reference group exhibits negligible sensitivity to the rotational speed, implying scarce production of additional hydrodynamics within the reference group.At lower rotation speeds,     www.nature.com/scientificreports/such as 15 rpm, the friction coefficient can be curtailed by 4.5% using optimal texture parameters, while it experiences a 0.7% increase in the reference group.Although the impact on the friction coefficient is minimal in both cases, it reveals that the effectiveness of texture design in enhancing hydrodynamic lubrication is limited at very low rotation speeds.Consequently, both optimal and initially designed texture parameters yield similar outcomes, suggesting that the surface texture design of the helical pair could better serve as a secondary lubrication and wear particle accommodation mechanism rather than providing additional hydrodynamic lubrication.At higher rotation speeds, the optimal texture design significantly curtails the friction coefficient, which stands at 0.0439 at 120 rpm.Compared to an untextured helical pair, the friction coefficient is reduced by 48.2% in the optimal texture design, while the reference group only achieves a 5.4% reduction.This translates to a 45.2% improvement in friction coefficient post-optimization, indicating that the friction coefficient is highly sensitive to texture parameters.A well-structured surface texture facilitates the generation of additional hydrodynamic lubrication at relatively high speeds.
Beyond 30 rpm, the optimal texture parameters remain almost unchanged and could be chosen as the optimal design for the helical pair, considering the limited influence of surface texture on friction coefficient at extremely low speeds.The optimal texture parameters are as follows: The number of radial dimples Z r is 2, the number of circumferential dimples Z c is 6, the dimple depth h d is 4-5.5 um, the axial ratio of the ellipse r t is 3, and the orientation of the major axis φ is 0 • .The geometric characteristics of the optimal texture parameters can be summarized as follows: a relatively small number of dimples in the radial and circumferential direction, implying larger dimple diameter and spacing, an appropriate dimple depth, a relatively flat ellipse shape, and an elliptical dimple shape with its major axis parallel to the sliding direction.These characteristics greatly enhance the hydrodynamic effect and significantly reduce the friction coefficient.The geometric characteristics of the optimal results align well with the study conducted by Yousfi et al. 41 , corroborated by experimental results.This is also evidenced by Fig. 11a, where compared to the reference group's relatively large number of deep, circular dimples shown in Fig. 12a, the smaller number of shallow, elliptical dimples with a sliding direction parallel to the major axis creates a larger converging wedge, thus generating more additional hydrodynamic pressure.

Effects of multiple parameters on friction coefficient
Given that a rotational speed of 120 rpm is a typical operating condition for the HHRA, and the most notable friction reduction effect is achieved at this speed following texture optimization, all subsequent parameter analyses will be conducted at 120 rpm.Moreover, the final values for the correlation parameter θ k in the Kriging model, after iterations, are [0.0625,0.0210, 4.0000, 2.6918, 0.7430].
Understanding the impact of each design variable on the optimization target is crucial.To intuitively and quantitively estimate the influence of each parameter on the friction coefficient, a popular global sensitivity analysis method, Sobol's method 42 , is utilized based on the built surrogate model.Figure 13a depicts the results where the bar lengths represent the total order sensitivity index and the first order sensitivity index, signifying the level of importance in affecting the friction coefficient.While the first order sensitivity index quantitively illustrates the significance of each variable, the total order sensitivity index conveys not just the effect of each variable on the target but also the interrelations among all design variables.As depicted in Fig. 13a, the first order sensitivity indices for the optimum axial ratio of the ellipse r t , the number of radial dimples Z r , the number of circumferen- tial dimples Z c , the dimple depth h d , and the orientation of the major axis of the ellipse φ are 0.006, 0.012, 0.024, 0.808 and 0.028, respectively.Corresponding total order sensitivity indices are 0.017, 0.022, 0.097, 0.919, and 0.115.The markedly higher first order sensitivity indices of the dimple depth demonstrate its major contribution to the friction coefficient of the helical pair under these conditions.The negligible difference between the first order and total order sensitivity of h d suggests that, with constant area density, the effects of interactions between h d and other parameters on the friction coefficient are minimal and largely mutually independent, underscor- ing the pivotal role of dimple depth in lubricating film formation.With increased dimple depth, both the local additional hydrodynamic pressure and lubricating film thickness decrease, leading to enhanced micro-contact  of asperities between the contact surfaces and subsequently, a rise in friction coefficient, and vice versa.Apart from dimple depth, the first order sensitivity indices of other parameters are quite marginal.However, the total order sensitivity indices of these parameters significantly surpass their first order counterparts, particularly for Z c and φ , where total order sensitivity indices are approximately four times their first order sensitivity indices.This indicates that interactions among parameters, excluding dimple depth, have a greater influence on the friction coefficient than individual parameters.In essence, individually modulating parameters apart from dimple depth may not significantly alter the friction coefficient.However, adjusting these parameters in unison could considerably affect the friction coefficient.Therefore, it is crucial to judiciously match parameter relationships to positively impact friction reduction.
While the aforementioned sensitivity analysis can assist in understanding the degree to which each design variable of the texture structure impacts the friction coefficient, it does not reflect the relationship of the friction coefficient with variations in the design variables.In addition, the accuracy of the sensitivity analysis, which is based on a surrogate model, necessitates further verification.Therefore, this study employs the Spearman rankorder correlation 43 to further investigate the relationship between the design variables and the objective function, and simultaneously verify the accuracy of the sensitivity analysis.Figure 13b exhibits the correlation coefficients of each design variable.It is found that the friction coefficient correlates with Z r and h d , and negatively correlates with r t , Z c , and φ .The correlation coefficient value of h d is greater than 0.8, indicating a strong positive correlation between the friction coefficient and the dimple depth.Consequently, a reasonable reduction in the dimple depth can effectively lower the minimum pressure.The absolute values of the correlation coefficients of the remaining design variables do not exceed 0.5.Therefore, it can be inferred that there is a lack of tight connection between the friction coefficient and the design variables, excluding h d .This is consistent with the conclusion derived from the sensitivity analysis, further verifying its accuracy.The aforementioned correlation analysis and sensitivity analysis offer an effective method for exploring the internal relationships of the coupled parameters.
Figure 14 presents the contour plots of the friction coefficient, where each subplot reveals the influence of two design variables on the friction coefficient while the remaining variables are held at their median values, derived from the averages of the respective upper and lower boundaries.The contour plots related to h d exhibit www.nature.com/scientificreports/ the maximum friction coefficient difference of 0.0223, surpassing that of other design variables.This suggests that h d adjustments can significantly alter the friction coefficient, establishing its pronounced influence.Conversely, other design variables individually exert a lesser effect compared to h d , particularly for Z c and r t .This pattern aligns with the results in Fig. 13a.Moreover, in contour plots associated with h d , the minimum friction coefficient value consistently appears near the first quintile of the h d range, approximately 4 mm.These parameters align with the optimal results from 5.2.The large gradient around this value shows that minor h d adjustments in this region can noticeably change the friction coefficient.In contrast, within the posterior two-thirds of the range, the small gradient signifies a lack of sensitivity of the friction coefficient to fluctuations in h d .This indicates that with other design variables fixed, the friction coefficient markedly decreases to the trough and then slowly ascends with an increase in h d .Additionally, in the φ − Z c contour plot, the maximum friction coefficient dif- ference is 0.0160, whereas the remaining plots unrelated to h d do not exceed 0.008.This suggests the significant influence of the combined effect of φ and Z c on the friction coefficient.The variables φ and Z c directly determine the circumferential distance between the dimples, significantly impacting the hydrodynamic effects.Increasing Z c and adjusting φ to favor the sliding direction reduces this distance, causing the additional hydrodynamic load capacity generated by each adjacent dimple to interact, thereby affecting the overall friction coefficient.In contrast, opposite adjustments to Z c and φ will increase the distance, reducing the interaction impact, which is directly reflected in the friction coefficient change.This principle is consistent with the results in Fig. 13a, where the total order sensitivity indices of φ and Z c are approximately four times their first order sensitivity indices.
From the above analysis, it is apparent that the friction coefficient is mainly influenced by the independent effect of h d , as well as the combined effect of φ and Z c .To further explore whether the combined effect of φ and Z c on the friction coefficient changes with variations in h d , other design variables are set to intermediate values.Figure 15 respectively illustrate the trend of friction coefficient changes with φ and Z c at different dimple depths h d .It can be observed that all surfaces exhibit a "saddle shape", i.e., they bulge at the sides and middle in the Z c direction and are concave at the sides and bulge in the middle along the φ direction.The consistency across all surfaces also suggests that the combined effect of φ and Z c is minimally influenced by h d .Furthermore, within all surfaces, the friction coefficient generally escalates with the increase of h d , which is in agreement with the previous analysis.

Conclusion
This paper introduces a modified efficient global optimization method predicated on the Kriging model, devised to optimize the surface texture of the helical pair in HHRA.Initially, a surrogate model was constructed by integrating the lubrication model of the helical pair with MEGO.Subsequently, to affirm the efficacy of the proposed MEGO in surface texture design, a comparative analysis involving MEGO and multiple prevalent global optimization algorithms was undertaken.Ultimately, surrogate-based optimization and analysis were executed with the goal of minimizing the friction coefficient.Based on the analyses and results, the ensuing conclusions are made: (1) Upon comparing various evaluation indices of the global optimization algorithms, the MEGO, in contrast to the EGO, Kriging + GA, BP + GA, and RBF +GA, demonstrates admirable proficiency in predicting the friction coefficient and optimizing texture parameters.This suggests that the MEGO can procure the optimal solution more rapidly while sustaining prediction accuracy.
(2) The dimple depth exerts the most substantial influence on the friction coefficient of the textured helical pair when the area density remains unchanged.Furthermore, the design parameters, excluding the dimple depth, interact mutually, exerting a more pronounced effect on the friction coefficient than any individual parameter, particularly in terms of the interaction between the orientation of the major axes of elliptical dimples and the number of circumferential dimples.
(3) Optimal texture parameters imply that large-sized, relatively flat elliptical texture features, with an appropriate dimple depth and the major axis oriented along the sliding direction, result in the most substantial reduction of the friction coefficient in the mixed lubrication regime.
(4) With the optimal texture parameters, the friction coefficient can be curtailed by as much as 45.2% compared with the initial design, underscoring the effectiveness of the proposed method in the surface texture optimization design process for superior tribological performance.

Figure 2 .
Figure 2. Variation of w iter with i iter .

Figure 3 .
Figure 3. Equivalent section of the interface of helical pair.

Figure 4 .
Figure 4. Cross-section of the screw tooth.

Figure 5 .
Figure 5. Equivalent section of the interface of helical pair with elliptical-shape dimples.

Figure 6 .
Figure 6.Flow chart of optimization design for the surface texture of helical pair.

Figure 11 .
Figure 11.Optimum design of surface texture.(a) Pressure distributions.(b) Top view of texture profile.

Figure 12 .
Figure 12.Initial design of surface texture.(a) Pressure distributions.(b) Top view of texture profile.

Figure 14 .
Figure 14.Multiple plots of friction coefficient with respect to design variables.

Figure 15 .
Figure 15.Relationship of µ cof with variations in Z c and φ at different h d . (a) h d = 5 um (b) h d = 10 um (c) h d = 15 um (d) h d = 20 um.

Table 1 .
Geometric and operating parameters.

Table 3 .
Computation time for several different grid densities.

Table 4 .
Results of the different global optimization algorithm for texture parameters.sample points is set to 200, the Kriging model at this juncture meets the accuracy requirements.As depicted in Fig. 9b, no changes are observed in the minimum friction coefficient values as the number of iteration steps increases when the number of sample points exceeds 130. From the previous analysis of the MEI, it is recognized that the MEI is capable of balancing local development and global exploration.No further reductions in the friction coefficient are observed when the number of sample points increases from 130 to 200.Consequently, it can be inferred that the proposed algorithm has potentially located the global minimum of the friction coefficient.

Table 5 .
Final optimization results at different rotational speeds.